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Abstract 

This paper presents a set of realizable second order models for boundary free turbulent 
flows. The constraints on second order models based on the realizability principle arc re- 
examined. The rapid terms in the pressure correlations for both the Reynolds stress and 
the passive scalar flux equations are constructed to exactly satisfy the joint realizability. 
All other model terms (return- to-isotropy, third moments and terms in the dissipation 
equations) already satisfy realizability (Lumley 1978, Shih and Lumley 19SG). To correct 
the spreading rate of the axisymmetric jet, an extra term is added to the dissipation 
equation which accounts for the effect of mean vortex stretching on dissipation. The test 
flows used in this study are the mixing shear layer, plane jet, axisymmetric jet and plane 
wake. The numerical solutions show that the new unified model equations (with unchanged 
model constants) predict all these flows reasonably as the results compare well with the 
measurements. We expect that these model equations would be suitable for more complex 
and critical flows. 

I. Introduction 

The second order closure scheme has been studied for over two decades now and it is 
playing an increasingly important role in the computation of turbulent flows, for example, 
in atmospheric turbulence and turbulent combustion. So far this method has achieved 
success in predicting many different flows (Launder et al^, Lumley et al^, Shih and 

Lumleyl I. * 3 * !, ShihW, Chen! 5 ] and Ettestad^ ). However, as Schumann^ and Lumley^ have 
pointed out, many of the second order model equations are not realizable because the 
models for the pressure correlations in the second moment equations do not satisfy the 
realizability constraints. They also pointed out that this may cause severe numerical dif- 
ficulties and produce unphysical results during a numerical computation; such as giving 

negative turbulent energy or producing correlation coefficients larger than unity in cer- 

tain critical situations. Hence, there is a strong need for improving second order model 
equations for both theoretical reasons and practical needs. In this paper we follow Shih 


1 


and Lumleyf 9 ! to directly impose the realizability principle (including joint realizability 
between the velocity and passive scalar) in the model development and obtain a set of 
model forms for the rapid part of the pressure correlations. In addition, an extra term is 
added to the dissipation equation which reflects the effect of the mean vortex stretching 
on the dissipation. 

In order to discuss various model terms appearing in the second moment equations, 
let us write down the exact equations for the mean quantities and various second moments 
for incompressible flows (including passive scalar): 


Ui,i = 0 

D <t Ui = ~~ p ,i ~ + vUi,jj 


D,F = -(fu j ) tj + 1 F t 


JJ 


D,tU{Uj = - [iiiUjU k + ~(p ( 2 hiiSj k + p^UjS'k)} 
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UiUkUj'k - UjU k U itk - Hp^Uj + Pj’ui) 

1- 


, 0 ). 


+ -(pWuij +pWuj'i ) - 2vu tk u ]k 
1- 


D't f ui = ~ [fuiUj + -pWfSij] 
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- UiUjFj - fujUij - -p { } ] f 
1- 


+ ~ 0' + 7 )/,>«.-,> 
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(4) 


(5) 
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where Ui,F and P are the mean velocity, mean scalar and pressure, and / and p arc their 
corresponding fluctuating quantities. Here, we have employed the summation convention 
on the indices and the following notations: ( ),, = D,( ) = Jy + U k ( ) ik . For the 

pressure fluctuation, the following Poisson equation must be satisfied: 




"f" ^ i-,3 u j,i U i,j U j,i 


(7) 


Based on the linearity of p, we have split the pressure fluctuation into two parts: p 1 - 1 and 
p( 2 \ called the rapid and slow pressures respectively: 


■ = 9 it- u ■ 


1 ( 2 ) 

-P,jj = 


(S) 

(9) 
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Note that the pressure gradient correlation terms involving p ^ (slow pressure) in the Eqs. 
(4) and (5) have been seperated into a pressure transport and a pressure strain (or pressure 
scalar gradient) correlation. 

For homogeneous turbulence, using the solution of the Poisson equation (8), we may 
write 


1 (i) 

-P y j u t = 
P 
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( 10 ) 

( 11 ) 


Equations (10) and (11) may give us a hint on how to construct models for the rapid 

pressure terms. 

For convenience, one often defines e = ( f — 7 They represent the 

mechanical and scalar dissipation rates respectively and must be modeled. Their transport 
equations may be written as follows 


D, t e = ^ 

- 2 


( 12 ) 


(13) 


Here, 'P and P j contain all the constructive and destructive terms of the dissipations (see 
section II), and q 2 = ujul. 

Lvimleyt 8 ! suggested that the slow pressure strain (or pressure scalar gradient) cor- 
relation term and the viscous dissipation term in the second moment equations can be 
combined and modeled together, because they are both related to only purly turbulent 
quantities. Therefore, we write 


-$ij€ = ~p (2 K u i,j + u J,i) ~ 2’' u i,kUj,k + 

P 6 


1 


-<Pi— = -p (2) /,, - (v + 7 


(14) 

(15) 


For the rapid terms, Eqs. (10) and (11), we need to model two integrals: 


Ipjqi ~ 4t r 
Iijk = ~ 4 ^ 



(1G) 

(17) 


In the next section, we will discuss the realizability principle and use it to construct 
the models for various unknown correlations in the second moment equations. In section 
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Ill, we test these models in the boundary-free shear flows: mixing shear layer, plane jet 
and axisymmetric jet, and plane wake. The calculations show that the present models 
(with unchanged model constants) predict all these flows reasonably as the results show 
good agreement with the measurements. 

II Second Order Closure 


Realizability Principle 

The concept of realizability was first introduced by Schumann^ 7 ! and Lumleyt 8 ! . The 
basic idea is that any non-negative turbulent quantities (say, turbulent energy compo- 
nents, intensity of scalar fluctuations, etc.) must remain positive during the evolution 
of turbulence, and Schwarz’ inequality between any turbulence quantities (say, between 
fluctuating velocities and scalars) must be satisfied at all the time. The exact turbulent 
equations derived from the Navier-Stokes equation, for example, Reynolds stress and scalar 
flux equations, possess these physical and mathematical properties, i.e. the solution of the 
exact turbulent equations satisfies realizability. However, modeled turbulence equations, 
obtained by various approximations from the Navier-Stokes equation, often violate real- 
izability and produce unphysical results. In fact, so far none of the second order closure 
models, with the exceptions of Shih and Lumleyt 9 ! and Ristorcellib 0 !, satisfy complete re- 
alizability (which includes joint realizability between velocity and scalar). In this section, 
we will discuss turbulence models based on realizability. To do that, we define a corre- 
lation tensor Dij — f 2 ujuj — f u, fuj which consists of the Reynolds stress and scalar 
flux. Lumleyl 8 !, Shih and Lumleyl 9 ! argued that the realizability principle stated above is 
equivalent to that of non-negativity of the eigenvalues of both the Reynolds stress tensor 
Rij = u t u j and the correlation tensor D l} . That is, those eigenvalues must remain posi- 
tive during the evolution of turbulence. The simplest way to ensure this realizability is to 
require that the first derivative of the eigenvalues should vanish and the second derivative 
should remain positive if the eigenvalues vanish, see Figure 1. If we designate the eigenval- 
ues of the Reynolds stress tensor and the correlation tensor by R aa and D aa respectively 
(no summation convention for Greek indices) we may write these realizability conditions 
as 


Raa 

-» 0 

if Raa 

-> 0 

(18) 

D,tD aa 

-* 0 

if D(x Q 

-» 0 

(19) 

D ^ttRaa 

IV 

o 

if Raa 

-► 0 

(20) 

D ,ttD aa 

IV 

o 

if D aa 

-> 0 

(21) 


Eqs.(18) and (19) are the necessary conditions for realizability, and Eqs.(20) and (21) 
together with (18) and (19) will provide the necessary and sufficient conditions for realiz- 
ability. For more details, see Lumleyf 11 ! and Shih et atf 12 l. 

To impose the realizability conditions on various model terms in the equations for 
the second moments, we need the equations for the eigenvalues of Rij and Dij. In other 
words, we need the equations for R,j and D tJ in the principle axes of R t) and Dij. In the 
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principle axes of i?,j , the Reynolds stress equation (4) becomes 


D it u 2 a = ~[u 2 a Uk + 2-p^u a S ak ] k - 2 ulU a ,a 

“t ^Up qlp 0 q a (^orar 4” 


( 22 ) 


where u ^ are the eigenvalues of Rij . Now, if we impose the realizability condition (IS) 
on Eq.(22), we may obtain a set of constraints for the model terms in the Reynolds stress 
equations: 

U p ,ql paqa 0 if < - 0 (23) 

(*.« + \)-e - 0 if u* -> 0 (24) 

(u 2 ut 4- 2-p (2 )u 0 i 0 fc) k -* o if u 2 — > 0 (25) 

Similarly, we may write an equation for Dij in its principle axes and impose the real- 
izability condition (19) to obtain the following constraints on the model terms appearing 
in both the Reynolds stress and the scalar flux equations: 


&p,q{f^fpaiqa f^a^apq) 

if D aci 0 


(26) 


2^fu a $ a -f 2 ($ QQ + ht-Ztf< - 0 

q z 6 

if D aa - 0 (27) 


2fu Q (fu a U k + -p {2) fS a k),k 
P 

- (f 2 u 2 u k + -p (2) u a s 0 k ) ik - v 2 (f 2 Uk),k -» o 
P 

if D aa - 0 (28) 

The constraints (23-28) on each model term in the Reynolds stress and the scalar flux 
equations will ensure that the model equations satisfy realizability. The models proposed 
by Lumleyf 8 ] for the third moments and the slow term 4>,j, and the model proposed by 
Shih and Lum- ley^ already satisfy the abovementioned constraints. What we need here 
are the models for the rapid terms: I p j q i and f jk . These terms are usually very important 
terms in the Reynolds stress and flux equations. Unfortunately, many existing models do 
not satisfy the conditions (23-28), and therefore, may produce unphysical results. 

Models of the Rapid Terms 

In the past it has been customary to express I p j q i as a simple linear function of h l} 
and I{j k as a simple function of fu t . We find that it is impossible for these forms to 
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satisfy realizability. The most general forms were first proposed by Shih and Lumleyf 9 ! 
(also see Shih et atf 12 l). Here, we adopt simpler forms for Ipj qi and I tJ k which are capable 
of satisfying realizability: 


l pjqi 


— fiqi^pj T &2 (.bpqbij T ^qj^pi) 


T Ot^Sqjbpj T O' 4 8 p j 6 


pj qi 


~{~ pqb{j ~f“ &ijbpq + ^jq^pi “h b pibqj^ 

~f~ &6^qibpj &7&pjbqi 

+ or s(6pqtfj + fiijbpg + &jqb 2 pi + bpibqj) 

4 otqbqibpj + o/io(b ))(l bjj + b q jbpi ) 

4 Giibqibpj + Oi\2bpjb 2 qi 

+ ^13 (b pg b 2 } + b t jb pq + b q jb pi + b pt b 2 j) 


(29) 


likj = (3iS ik 0uj + /? 2 ( SijOuk 4 S jk 6ui) + fa b tf; duj 
4 P*{ bijdui . + bjkBiti) + &(M* P + 8kjbip)0u v 
T Pq k bjp Oxip T fiqbifcbjpQup 
4 Ps (bijbkp + bkjbip)6u p + P$b 2 ik 9uj 
4 - Pio(b 2 jOu k + b 2 k 9u,j ) + Pi\9ikb 2 p9up 
4 P 12 (bij b kp 4 b j k blp)9up T P\^b ik bjp9xi p 

4 P\a bikb'jpOiip + /3is(b 2 jb kp + b 2 k jbjp)9 ti p (30) 

where bij = is called the anisotropy tensor of the Reynolds stress. The coeffi- 

cients and Pi in Eqs. (29) and (30) are, in general, functions of the invariants of b tJ and 
Dij. However, for passive scalar turbulence, cv; should be only a function of the invariants 
of bij. These coefficients need to be determined. To achieve this, we recall the definitions 
of Ipjqi and Iijki he., Eqs. (16) and (17), and find that they have the following properties: 

Ipjqi = I jpqii Ipjqi ~ Ipjiqilijk — I jik (31) 

Ippqi = UqUi, Ipkqk — O^Ippk — f Kki likk = 0 (32) 

We notice that Eqs. (29) and (30) already satisfy the condition imposed by Eq. (31). If we 
use the condition (32) and the realizability constraints (23) and (26), we may determine 
the limiting values of all the coefficients. The final expressions are surprisingly simple: 

- M,i) 

^( Iqibpj bpjbqi) T O ] ( 8pq bij 4 8ijb pq f 8 q jb p i 
11 4 

4 Opibqj T^bqiOpj — ~Opjbqi) 
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(33) 


"4 6?2 ( 2 ^ pj b q j 3 bpgbjj 9bqjbp{ T b q I bpj ) 


lijk — ^bijOuf; 4" bjl c Oui') 

+ CoibijOUk 


+ C D2{bik@ u j + bjkOui) + CD^bijbklOui, 
where the limiting values of the coefficients are 

1 _1_ 

10 ’ 


- “To’ 02 


Cd, = h Co, = C m = i 


(34) 


(35) 


The last line in (33) and the last two lines in (34) represent the non-linear contribution 
and if neglected, the linear models used by various other workers will be recovered. It is 
important to note that the above values of the coefficients oj, 02, C o\ • C ir> and Co 3 
are their limiting values at the realizability limit, i.e. when u a u Q , D aa — * 0. For general 
turbulent flows it a it a and D oa are not zero and hence the values of the coefficients may 
deviate from their limiting values. They are, in general, func tions of the invariants II and 
III for a\ and 02, and the invariants formed from b tJ and /it, for Coi- Some guidance can 
be obtained by inspecting the following two useful parameters (see Lumley ^ and Shih 
and Lumley ^): 


F = 1 + 27 III + 911 

_ „ ,, 27 l2 9 

Fd — 9 d tJ - -d- + — , 

where 

= m = hi 

, _ P UjU ~ - /it, fUj 
13 P u~TTi - /it, /it ,■ 

It can be shown that both F and Fd lie between 0 and 1, and particularly 


(30) 

(37) 


F — * 0 when u a u a — > 0 
Fd — * 0 when D aa — > 0 

By using this information, it is convenient to write 


= -^ (1 + Apa) 

(38) 

<■2 = ^(1 +BF°) 

(39) 

C m =h + Ci Fi 

(40) 

c D 2 = ^ + c,f; 

(41) 

Cm = \ + C,FS 
0 

(42) 
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where A, B, Cj , C 2 and C 3 are adjustable constants but a is not as arbitrary as it might 
seem at first look. The conditions (20) and ( 21 ) suggest a = £ (see Lumleyf 1 ^ or Shih 
and Lumley f 3 J). In the limiting case these coefficients reach their limiting values. Shih et 
al. took A = 0.8 and B = 0.0 which fit the DNS data quite well. Shih et al. f 14 k 15 l set 
C\ — Ci — C 3 = 0 in their computations. 


Models of Other Terms 

The focus of this paper is on the calculation of the velocity field in the boundary- 
free flows. Therefore, here we list only the related models for the pressure transport, the 
third moments, and the return-to-isotropy terms in the second moment equations and the 
models for the dissipation equation. All these models were proposed by LumleyW , and 
they satisfy the realizability conditions discussed in the section II. 


Pressure transport term: 


-pl 2 l U{ = Cq 2 Ui 


where C is a constant, and Lumleyf 8 ! suggested C = 0 . 2 . 

Third moments: 

_J_ r. ... , 

UiUjUk — 

+ Uj^(uiUk),p + u'iU p (ujUk),p] 

(3 - 2 

+ 1 Sijq 2 Uk + S lk q 2 Uj + 8 jk q 2 u t \ 


q 2 u k 


3 q 2 


4/3 + 10 ? 


— [u k u,,q 2 p + 2 u v u q u k u q ] 


(43) 


(44) 

(45) 


Return-to-isotropy term: 

= (3b ij 


(3 = 2 + exp( - + 80. 1 In [l 

+ 62.4 (— II + 2.3 III )} } (i + 3/17 + //) 


Model terms in the dissipation equation: 


eu k = - 


V/7 _ 

5(4/? + 10) 


r UfcUg Hq Up -| 

,p[u k u p + 2 % J 


$ = 0o + 01 ^bijUij 

0° = y + 0.98[exp(J^)][l 
-0.331n(l -55 //)] + 0 cor 


(46) 

(47) 

(48) 

(49) 

(50) 
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where ip\ — 2.4 is a model constant. The term j/> cor is similar to that proposed by Popef 16 !, 
which represents the effect of mean vortex stretching on the dissipation: 

Vw = 1.25(1 - F)° ‘ - UiMVi.k 

- + U, it ) (51) 

For isotropic turbulence (F = 1), xp cor becomes zero. For planar flows, tf> cor is also 
zero because there is no mean vortex stretching. With this extra term in the dissipation 
rate equation, the spreading rates are predicted very well for both the planar and the 
round jets. 


Ill Boundary- Free Shear flows 

This section present s the results of calculations for some boundary free turbulent shear 
flows including the two dimensional mixing shear layer, planar jet, axisymmetric jet, and 
two-dimensional wake. The numerical solutions were obtained by simultaneously solving 
the set of equations for the mean momentum ?7j, Reynolds stress TTpUj and dissipation 
e. All the model constants in the equations remain the same for all the test flows. For 
these thin shear layer flows, we have adopted the boundary layer approximations, and, 
therefore, the modeled equations are parabolic (we have kept the viscous diffusion terms 
in the modeled equations). The numerical scheme is based on the method of Patankar and 
Spaldingl 17 ^ 18 !. 

The boundary conditions imposed on these flows arc the following: for the mixing 
layer, the upper and lower free stream velocities have specified values, the derivatives with 
respect to the transverse direction of all other variables at the upper and lower boundaries 
have been set to zero. For the jets, the free stream velocity is zero, and the turbulent 
shear stress is set to zero at the center line of the jets. The transverse derivatives of all 
other variables at the boundaries are zero (including the mean velocity at the center line). 
For the two-dimensional wake, the free stream velocity has a specified value. All the other 
boundary conditions are the same as for the jets. 

The initial profiles of all the quantities are arbitrary smooth profiles. The calculations 
show that all the solutions had reached self-preservation. All the figures presented here 
are from the far field solutions. 

Figures 2, 3 and 4 show the profiles of the mean velocity, turbulent shear stress and 
energy components for the mixing layer. The experimental data were taken from Bradshaw 
et alt 19 !, Castro^ 20 ] and Gutmark & Wygnanskil 21 !. The computed mean velocity is in 
very good agreement with the measurements. The shear stress profile is also satisfactory. 
The experimental data of the energy components possess scatter but the model shows 
reasonable agreement with the experiments. The spreading rate (defined as dh/dx , h 
being the lateral distance between the positions where the velocity is 90% and 10% of the 
free stream) is calculated to be 0.13, which is also within the experimental scatter. 

Figures 5, 6, 7, 8 and 9 show the profiles of the mean velocity, shear stress and en- 
ergy components for the planar jet. The measurements were taken from Bradburyf 22 ^, 
Heskestadt 23 ! , and Gutmark & Wygnanskif 21 !. The prediction of the mean velocity is in 
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good agreement with the measurements except in the region ^ y-A r 0 ) > ^-1, where the 
model gives a slight underestimation. The calculations of the shear stress, streamwise and 
transverse energy components are within the scatter of the experimental data. However, 
the lateral energy component is apparently overestimated with respect to _the measure- 
ments. It should be pointed out that this set of measurements shows w 2 < v 2 and this is 
not consistent with the measurements for other shear flows (say, the mixing layer, wake 
and axisymmetric jet). The calculated spreading rate (defined by dY 5 /dx , Y 5 being the 
position where the velocity is the 50% of the centerline velocity of the jet) is 0.11 which is 
very close to the measurements. 

Figures 10, 11, 12, 13 and 14 show the profiles for the mean velocity, shear stress 
and energy components in the axisymmetric jet. The calculations are compared with the 
measurements of Abbiss et alt 24 !, Wygnanski Sc Fiedler^ 25 ! and Rodi^ 26 *. The predictions for 
all the quantities are in good agreement with the measurements. The calculated spreading 
rate (defined as the same as in the planar jet) is 0.09 which is also very close to the 
measurement s. 

Finally, figures 15, 16 and 17 show the profiles for the mean velocity, shear stress and 
energy components in the two-dimensional wake. Usually the 2D wake is a strongly non- 
equilibrium flow. In our calculation, it takes more time for the solution to approach self 
preservation as compared with the solutions of the mixing layer and jets. The predictions 
of various quantities agree reasonably well with the measurements at the far field of the 
wake, even though the measured wakes are probably not becoming self-similar yet. 

From the above calculations and comparisons, we conclude that the model based 
on the realizability concept performs quite well for typical boundary-free turbulent shear 
flows. The modeled equations are realizable and will not produce unphysical results, and 
therefore, we expect that the present model would be suitable for more complex flows. 
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Fig. 1 Realizability condition when eigenvalues vanish. 






Fig. 2 Mean velocity profile in 2-D mixing layer. 

(lie free stream velocity, Y 5 : the position where U = 
hUmax- O : Bradshaw, et atf 19 !, — : Present model. 



Fig. 4 Normal stress profiles in 2-D mixing layer. O. A, V 
Castrd 20 !, Q.A.V: Gutmark k Wygnanskit 211 , The 
lines represent the present model. 
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Fig. 10 Mean velocity U/Us profile in axisymmetric jet. 
Us' the centerline mean velocity. 0 : Abbiss et al^, 
X: Wygnanski k FiedleJ 25 ^ — : Present model. 



Fig. 11 Shear stress uv/U$ profile in axisymmetric jet. 
Q: Rodi^ 26 ^ X : Wygnanski k Fiedler^ 25 !, A: Abbiss et 
al, — : Present model. 



Fig. 12 Normal stress u 1 jU\ profile in axisymmetric 
jet. Legend as in Fig.ll. 
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Fig. 13 Normal stress v 1 jU\ profile in axisymmetric 
jet. Legend as in Fig.ll. 



Fig. 14 Normal stress w 1 /U\ profile in axisymmetric 
jet. Legend as in Fig.ll. 
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Fig. 15 Mean velocity profile (U\ — U)/Us in 2 -D wake. 
U\: the free stream velocity, Uo’ the centerline mean 
velocity, Us = U i — Uo, Yg: the position where U\ — 
U — O. 6 O 5 . X: Chevray &; Kovasgnay^ 27 !, — : Present 
model. 


Fig. 17 Normal stress profile in 2-D wake. Chevray & 
Kovasynayl 27 !; 0 : A: y2 /^s^ V ; w2 /^s • The 

lines represent the present model. 
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Fig. 16 Shear stress uv/Ug profile in 2-D wake. Legend 
as in Fig. 15. 
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